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ABSTRACT 

We  present  a simple  matrix  representation  of  the  Berlekamp- Massey  algorithm  for  the  minimum  partial  re- 
alizations problem,  and  show  how  pivoting  can  be  added  to  the  algorithm  to  improve  numerical  accuracy  of  the 
method. 


1.  Introduction 


The  problem  of  minimal  realization  of  linear  dynamical  systems  from  input/output  data  has  much  practical 
importance.  Many  of  realization  procedures  that  have  been  developed  rely  on  the  solution  of  a Hankel  system 
of  linear  equations.  The  Berlekamp- Massey  (BM)  algorithm  [1],  [6]  is  a fast  Hankel  linear  system  solver  which 
originated  in  the  field  of  coding  theory.  The  algorithm  is  little  known  in  the  scientific  computing  community.  One 
reason  for  its  obscurity  may  be  that  the  algorithm  seems  to  lack  a natural  representation  in  matrix  forms.  Attempts 
to  alleviate  this  situation  can  be  found  in  Rung  [4]  and  Jonckheere  and  Ma  [3].  In  this  paper,  we  give  a related  but 
perhaps  simpler,  way  to  present  the  algorithm.  We  show  how  our  presentation  leads  to  a pivoting  strategy  that 
improves  the  numerical  accuracy  of  the  computation.  What  is  more,  unlike  other  pivoting  schemes  for  Hankel  and 
Toeplitz  matrices,  our  new  algorithm  never  requires  more  than  0(n 2)  operations. 

Throughout  this  paper,  unless  otherwise  stated,  all  matrices  are  nxn  and  all  vectors  have  n elements.  Wherever, 
convenient,  we  will  use  upper  case  Latin  letters  to  denote  matrices,  lower  case  Latin  letters  to  denote  vectors,  and 
lower  case  Greek  letters  to  denote  scalars.  A Hankel  matrix  H has  the  form: 
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and  we  are  interested  in  solving  the  matrix  equation: 

Hx  = b. 

By  rearranging  its  columns  or  rows,  the  Hankel  system  can  be  transformed  into  a Toeplitz  system.  For  example, 
we  can  re-order  the  columns  of  H from  the  last  to  the  first,  and  get  the  Toeplitz  system  of  equations: 
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Similarly,  we  define  the  Yule- Walker  problem  for  the  Hankel  matrix  to  be: 


(1.4) 


Many  algorithms  for  solving  (1.3)  have  been  proposed,  but  most  of  them,  e.g.,  Levinson  [5],  may  fail  to  calculate 
an  accurate  solution  if  the  Toeplitz  matrix  has  ill-conditioned  principal  submatrices.  Interestingly,  our  numerical 
experiments  indicate  that  our  new  algorithm  may  still  work  very  well  under  these  circumstances.  There  is  much 
recent  interest  to  introduce  pivoting  to  Toeplitz  algorithms;  see,  e.g.,  [2]. 

This  paper  is  organized  as  follows.  Section  2 describes  how  one  solves  a Hankel  matrix  equation  via  the  BM 
algorithm.  Section  3 explains  how  the  BM  algorithm  triangularizes  a Hankel  matrix  that  is  strongly  nonsingular. 
Section  4 presents  our  new  numerical  pivoting  strategy.  Section  5 considers  the  case  of  a general  Hankel  matrix. 
The  last  three  sections  contain  examples  that  detail  our  numerical  experience. 
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2.  Solving  a Hankel  Matrix  Equation 


In  Section  3,  we  will  show  how  the  Berlekamp-Massey  algorithm  constructs  an  upper  triangular  matrix  72  to 
reduce  a Hankel  matrix  H to  a lower  triangular  matrix  L : 


HR  = L. 


(2.1) 


The  triangular  matrix  R also  has  a unit  diagonal.  From  this  factorization,  the  Hankel  system  (1.2)  can  be  easily 
solved.  Now, 


(. HR)t  = Lt  and  HT  = H 


imply  that 


RtH  = Lt. 


Multiplying  both  sides  of  (1.2)  by  72T,  we  get 

Ltx  = RTb.  (2.2) 

Hence  we  first  apply  RT  to  6,  and  then  solve  the  triangular  system  (2.2).  So,  if  the  factorization  (2.1)  is  available, 
a total  of  n2  multiplications  is  required  to  solve  (1.2).  It  is  worthwhile  to  point  out  here  that  the  matrix  R needs 
not  be  upper  triangular.  Even  if  R were  a dense  matrix,  one  could  still  solve  (1.2)  via  (2.2),  albeit  at  a cost  of  an 
additional  n2/ 2 multiplications.  When  we  introduce  pivoting  in  Section  4,  we  may  destroy  the  triangularity  of  72. 


3.  Hankel  Matrix  Triangularizafcion 


For  this  section  and  the  next,  we  assume  that  the  Hankel  matrix  H is  strongly  nonsingular,  i.e.,  that  all  its 
principal  minors  are  nonzero.  This  assumption  simplifies  our  presentation,  and  will  be  removed  in  Section  5.  For 
convenience,  we  need  a “shift-down”  matrix: 
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where  Jn_i  is  the  identity  matrix  of  order  (n  - 1).  Thus, 
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We 

proceed 


now  show  how  the  BM  algorithm  computes  columns  of  the  two  matrices  72  and  L of  (2.1)  recursively.  We 
by  induction,  and  use  the  usual  notation  representing  the  columns  of  the  three  matrices: 


H — {h\y  /i2>  ■ ■ ■ > hn), 
72  = (ri,r2,---,rn), 
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The  first  two  columns  of  the  matrices  R and  L are  readily  available: 


i\\ 


r2  = 


/ -mhi  \ 


1 

o 


w 


V 0 ) 


(3.2) 


and 

h = hi  , l2  = h2  - (W^i)  (3.3) 

Hence  the  top  element  of  l2  equals  zero.  Suppose  now  that  the  four  columns  r;,  l:  and  h+ 1 have  already  been 
computed,  and  that 

H ( ri  ri+i ) = ( h O’-H ) * (3-4) 


Let  us  denote  the  elements  of  r;*+1  and  h+i  by 
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From  the  assumption  of  strong  nonsingularity  we  are  assuming  that  ^ 0.  Also,  let 


n+ 2 = Zrj+U 

(3.5) 

and 

h+ 2 = ^+2- 

(3.6) 

That  is, 

H ( r;  ri+i  f/+2  ) = ( (/  O'+i  0+2 ) ■ 

The  new  vector  lj+2 

is  easy  to  compute.  From  (3.1),  it  is  seen  that 

(7+2  = %Tlj+ 1 +^en, 

(3.7) 

where 
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and  en  denotes  the  last  column  of  the  n x n identity  matrix.  In  words,  the  vector  /J+  2 is  formed  by  “upshifting” 
each  element  of  /;+i  by  one  slot  and  placing  the  scalar  £ in  the  n-th  position.  A picture  is  worth  a tho-isand  words! 
Hence  we  get 
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We  now  zero  out  the  two  leading  nonzero  elements  of  lj+ 2 by  using  appropriate  multiples  of  the  leading  elements 
of  lj  and  (*+1  • That  is,  we  post-multiply  the  n x 3 matrix  of  (3.9)  by  the  two  3x3  elimination  matrices  and 
E^\  where 

(l  0 m!0)\  (l  0 0 \ 

E[0)  =01  0 , 40)  = I 0 1 m(20)  ) , (3.10) 

Vo  o i / Vo  o i / 


and  n»|°\  mi0)  denote  the  corresponding  multipliers.  Finally, 

( h h+i  (j+2  lj+ 1 lj+2  )E[0)El°\  (3.11) 

and 

( Tj  r;-+ 1 rj+2  ) *—  (rj  rj+ 1 rj+2 ) E[0)E^0\  (3.12) 

Note  that  the  (j  + 2)-nd  component  of  rj+2  is  nonzero.  Hence  from  the  strong  nonsingularity  assumption,  the 
(j  + 2)-nd  component  of  lj+2  is  also  nonzero.  Thus  the  multipliers  are  well-defined. 

Let  us  perform  an  operation  count.  The  time-consuming  steps  include  the  calculation  of  the  inner  product  £ ( j 
multiplications  ),  and  the  multiplication  of  a scalar  into  the  four  vectors  rj , rJ+1,  and  lj+\  ( 2n  multiplications  ). 
Hence  a total  of  5n2/2  multiplications  is  required  to  compute  the  decomposition  (2.1). 

4.  Pivoting 

One  may  have  noted  that  the  magnitude  of  two  multipliers  mi  and  m2  of  (3.10)  can  be  arbitrarily  large.  In 
response,  we  propose  a simple  scheme  of  eliminating  the  leading  nonzero  element  of  either  l3  or  l3+ 2,  using  either 
4°>  or  23^,  respectively,  where 

/ 1 0 0\ 

4"  = 0 1 0 . (4.1) 

\mj1)  0 1/ 

The  important  point  is  that  either  or  must  be  at  most  one  in  absolute  value,  in  order  to  keep  the 
overall  process  stable.  We  thus  choose  either  E |°)  or  e[1^  to  achieve  a better  numerical  accuracy.  Our  approach  is 
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The  location  (i,  j)  of  the  multiplier  represents  the  non-zero  leading  entry  of  column  j is  to  be  eliminated  by  that 
of  column  i. 

The  column  updating  proceeds  essentially  as  before: 

(h  h+i  ';+2  )<-(/;■  h+ 1 h+2  )EiE2,  (4.3) 


(ri  rj+1  rj+2)*-(rj  rj+l  rj+2)ExE2)  (4.4) 

where  E\  equals  E^  or  23,1),  and  E2  equals  E E^\  E^  or  E^\  Two  important  observations  are  as  follows. 
First,  the  resultant  matrix  L stays  lower  triangular,  but  the  previously  upper  triangular  R may  have  gained  two 
nonzero  subdiagonals.  Second,  our  pivoting  scheme  increases  the  number  of  multiplications  by  only  0(n),  i.e  , the 
total  number  of  multiplications  is  still  5n2/2  + 0(n). 
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5.  General  case 


Recall  that  under  the  strong  nonsingularity  assumption,  we  knew  that  the  (j  4-  l)-st  element  of  0+i  must 
be  nonzero.  We  now  remove  the  assumption  that  the  matrix  H is  strongly  nonsingular.  During  the  elimination 
process,  we  may  get  additional  leading  zero  elements  in  l3 j i . For  our  discussion  in  this  section,  let  us  assume  that 
both  (j  4-  l)*st  and  ( j 4-  2)-nd  elements  are  zero  but  that  the  (j  4-  3)-rd  element  is  nonzero.  Hence  the  procedure 
described  in  Section  3 would  not  work  because  there  is  a gap  in  the  nonzero  structure  in  (3.7).  Now,  let 
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Define  some  new  vectors  by 
and 

That  is, 


rj+2  = Zi*j+ 1,  fj+  3 = Zrj+ 2,  0+4  = Z?j+  3, 

0+2  = ^0+2,  0+3  = Hrj+ 3,  0+4  = #0+4- 


H(rj  rj+1  fi+ 2 ?i+ 3 ri+4)  = (/i  ZJ+1  lj+2  lj+3  Zi+4). 

From  (3.1),  the  new  vectors  Zj+2,  O' +3  and  0+4  are  calculated  by 

0+2  = ZTlj+l  +^len. 


0+3  = ^T0+2  + fc^n, 

0+4  = ZT  lj+3  + ^3em 

where 

£1  = tfn+lPl  + Pn+2P2  H h Vn+j+lPj+l, 
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Indeed,  the  million  words  picture  looks  like: 
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As  described  in  Section  4,  we  would  like  to  pivot  and  eliminate  the  above  matrix  so  that  each  row  contains  an 
unique  pivot  element.  The  elimination  matrices  in  this  case  are  four  5x5  matrices,  each  with  all  l’s  on  the  diagonal 
and  a multiplier  in  the  (i,  j)  location.  Same  as  before,  the  jth  leading  non-zero  entry  is  to  be  eliminated  by  the  :th 
column,  and  all  the  multipliers  in  the  elimination  are  less  than  one  in  absolute  value. 

In  practice,  we  work  with  finite-precision  arithmetic,  so  an  exact  zero  would  hardly  happen.  In  order  to  tell 
if  we  are  getting  any  additional  zeros  in  the  column  of  /j+i,  we  need  to  choose  a threshold,  such  that  any  number 
smaller  (in  absolute  value)  than  the  threshold  is  regarded  as  a zero.  If  this  is  the  case,  we  will  then  apply  the 
technique  in  this  section  to  deal  with  the  situation. 

6.  Numerical  Examples 

We  consider  the  Hankel  matrix  equation  (1.2)  and  the  corresponding  Toeplitz  matrix  equation  (1.3).  We 
compare  three  procedures:  the  BM  algorithm  for  (1.2),  our  new  pivoted  BM  algorithm  for  (1.2),  and  the  Levinson 
algorithm  for  (1.3).  We  construct  two  sets  of  examples,  the  first  where  BM  would  fail  and  the  second  where 
Levinson  would  fail.  Specifically,  we  tinker  with  the  2x2  leading  submatrices  of  the  Hankel  and  the  corresponding 
Toeplitz  matrices: 


In  Example  1,  the  submatrix  H ^ is  ill-conditioned  but  the  submatrix  T ^ is  not,  while  in  Example  2,  the  situation 
is  reversed.  In  Example  3 both  submatrices  are  ill-conditioned.  Whereas  the  BM  algorithm  fails  in  Examples  1 
and  3,  and  the  Levinson  algorithm  fails  in  Examples  2 and  3,  our  new  algorithm  works  well  on  all  three  sets  of 
equations. 

For  all  examples  in  this  paper,  we  choose  the  left  hand  vector  b such  that  the  solution  vector  x = 
(1  1 •••  1)T.  To  compare  the  algorithms,  we  calculate  ||  x - x ||2,  where  x denotes  the  computed  solu- 
tion. We  ran  our  examples  using  MATLAB  on  a Sun  Sparc  station.  In  this  section  we  choose  as  a threshold, 
e ||  H H2 j where  e («  2.22  ■ 10”16)  denotes  the  machine  precision.  We  use  k(M)  to  denote  the  condition  number 
with  respect  to  the  2-norm  of  a matrix  M. 


Example  I.  This  example  shows  why  pivoting  is  necessary  for  the  Berlekamp-Massey  algorithm.  Let 
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and  H[2^  = 


(2) 

The  submatrix  of  H\  is  ill-conditioned  when  6 is  small,  and  is  singular  when  6 is  zero.  As  expected,  the  BM 
algorithm  delivers  worse  accuracy  as  we  decrease  the  size  of  5.  The  matrix  H\  is  well  conditioned,  with  k(H\)  = 5.6. 
However,  the  BM  algorithm  determines  an  L that  is  ill-conditioned,  contributing  to  the  loss  in  accuracy  when  one 
solves  (1.2)  via  (2.2).  On  the  other  hand,  our  new  algorithm  computes  a very  well-conditioned  L . 


Table  1.  Error  Behavior  for  Example  1 
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Example  3.  This  is  an  example  where  both  BM  and  Levinson  algorithms  fail  because  the  submatrices 
and  lj2)  are  ill-conditioned: 
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Again,  here  the  Hankel  matrix  is  well-conditioned,  with  *(#2)  « 31.  Our  new  algorithm  calculates  a “good”  L, 
but  the  BM  algorithm  determines  an  L that  is  ill-conditioned. 


Table  3.  Error  Behavior  for  Example  3 
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The  three  examples  have  shown  how  our  pivoting  scheme  works  better  than  the  other  two  conventional  al- 
gorithms, However,  when  the  size  of  the  matrix  H increases,  roundoff  errors  may  accumulate  so  that  it  becomes 
hard  to  define  a numerical  zero.  If  we  choose  a small  threshold,  such  as  the  one  we  have  used,  then  the  condition 
number  of  any  principal  submatrix  may  become  as  large  as  the  inverse  of  the  threshold  and  a significant  loss  in 
accuracy  may  occur.  From  our  experiments,  we  observed  that  the  accuracy  of  the  solution  is  proportional  to  the 
largest  condition  number  of  any  principal  submatrices. 

On  the  other  hand,  if  the  threshold  is  large,  we  effectively  work  in  a lower  precision,  so  the  factorization  will 
have  limited  numerical  accuracy.  Therefore,  in  the  next  section,  we  experiment  with  a compromised  threshould 
value.  To  compensate  for  the  loss  in  accuracy  due  to  this  choice  of  a larger  threshold,  we  adopt  iterative  refinement 
at  the  end. 
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7.  Further  Examples 


In  this  section  we  show  how  our  pivoting  scheme  works  when  the  size  of  the  matrix  H increases.  We  observe 
that  with  an  increase  in  dimensions  there  is  a danger  of  underflow.  Hence  some  form  of  normalization  is  required. 
Our  choice  is  to  normalize  L to  make  the  diagonal  elements  all  ones,  so  that  underflow  can  be  avoided. 

We  construct  our  matrices  from  the  Toeplitz  examples  in  Sweet’s  paper  [7],  and  select  an  iterative  refinement 
scheme  for  improving  the  accuracy  of  the  initial  solution 

1.  Compute  rM  = Hx W - b. 

2.  Solve  LTy  = RTr^\ 

3.  Update  = xW  — y . 

The  criterion  for  ending  the  iterative  refinement  is  when 

||  rW  ||2<  10  e- 1|  if  ||2 -||  *(’•>  ||2. 


The  threshold  for  the  examples  in  this  section  is  chosen  as  10  ■ y/e- 1|  H |b- 


Example  4.  We  pick  an  example  to  show  how  iterative  refinement  improves  our  solution,  and  how  the  number 
of  refinements  is  affected  by  the  conditioning  of  principal  submatrices.  The  order-6  Hankel  matrix  is 

1 195/14  + 5 8 \ 

i/14  + 5 8 4 

8 4 -34 

4 -34  5 ' 

-34  5 3 

5 3 1 / 

For  5 smaller  in  the  magnitude  than  0.5,  the  matrix  is  well  conditioned  with  k(/L|)  less  than  100.  The  threshold  is 
approximately  equal  to  6 • 10-6.  For  5 = 0 the  order-3  principal  submatrix  is  singular.  Starting  with  5 = 10“2  and 
then  decreasing  it , we  can  make  this  submatrix  progressively  worse  conditioned  without  significantly  changing  the 
condition  number  of  HA. 


H 4 = 


3 

2 

6 

1 

195/14  + 5 

8 


2 

6 

1 

195/14  + 5 
8 
4 


6 

1 

195/14  + 5 
8 
4 

-34 


Table  4.1.  Error  Behavior  for  Example  4 


i 


BM 

Pivoted  BM 

BM 

Pivoted  BM 

5 

10~2 

10"2 

10-4 

10-4 

1.36  • 10-3 

1.36- 10“3 

1.36  • 10"5 

1.36 -10"5 

<L) 

2.29  • 107 

197 

2.29  • 10u 

197 

k(R) 

2.29  • 107 

182 

2.29  • 10“ 

182 

||x-x«»  lb 

1.06  • 10"8 

1.33 -10“ 13 

6.59  • 10"5 

2.94  - 10-11 

Hz-xW  lb 
||*-x(2>  lb 

7.61  • 10-16 

7.84  ■ 10~10 
5.27  • 10-15 

7 

O 

fH 

CO 

CO 

# Re  fine. 

1 

0 

2 

1 
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Table  4.2 . Error  Behavior  for  Example  4 (Continued) 


BM 

Pivoted  BM 

BM 

Pivoted  BM 

8 

10~5 

10"5 

1(T6 

10-® 

<w43)) 

1.36  ■ Hr6 

1.36  • 10"6 

1.36  • 10"7 

1.36  • 10-7 

k(L) 

2.29  • 1013 

197 

214 

122 

k(R) 

2.29  • 1013 

182 

202 

53 

||z-x(0)  ||2 

1.23  • 10"2 

1.12  • 10-1° 

1.47  • 10"7 

2.11  • lO"7 

||*-xW  lb 

5.93  • 10"6 

3.84  • 10"15 

8.89  • 10-14 

1.55  • lO"13 

II  x - x<2)  ||2 

2.71  • 10"9 

1.37  • 10“15 

4.29  • 10“15 

II  X - z(3)  lb 

1.24  • 1(T12 

II  x - *(4)  ||2 

3.67  ■ 10“15 

#Refine. 

4 

1 

2 

2 

Table  4.3.  Error  Behavior  for  Example  4 (Continued) 


BM 

Pivoted  BM 

BM 

Pivoted  BM 

8 

k(L) 

k(R) 

||x-x(0)  ||2 

II  lb 

#/?e/ine. 

10"8 

1.36  • lO”9 
214 
202 

1.47  • 10"9 
2.51  • 10“15 
1 

10"8 

1.36  • 10"9 
122 
53 

2.11  • 10“9 
3.73  • 10"15 
1 

iO-io 

1.36 -10"11 
214 
202 

1.47  • 10“n 
5.87  • 10"16 
1 

10-10 

1.36  • 10_u 
122 
53 

2.11 -lO-11 
3.92  • lO"15 
1 

Note  that  the  pivoted  BM  algorithm  always  produces  factors  R and  L that  are  better  conditioned  than 
those  produced  by  the  BM  algorithm  without  pivoting.  As  a consequence,  the  first  approximation  to  the  solution 
computed  by  the  pivoted  BM  algorithm  is  more  accurate  than  that  computed  by  the  BM  algorithm.  When  the 
the  smallest  singular  value  of  the  order  3 principal  submatrix  becomes  smaller  than  the  threshold,  both  algorithms 
behave  in  a similar  way. 

Example  5.  We  construct  a 13  x 13  Hankel  matrix  H$  whose  first  row  is  given  by 

T)i-13  = ( -15, 10, 1,  -7,  -2,  -5,  -14.2766,  -25.5087,  -48.8789,  -96.8384,  - 188.8878,  -1,5), 


and  the  last  column  by 


7713-25  = ( 5, 1,  -3, 12.755,  -19.656, 28.361,  -7,-1, 2, 1,-6, 1,  -0.5  f , 
the  threshold  is  5.35  • 10“5. 

The  matrix  is  well  conditioned  in  that  k(Hs)  = 89.0,  but  it  contains  five  consecutive  ill-conditioned  submatrices 
to  H$S\  i.e.,  orders  4 through  8.  The  smallest  singular  values  of  these  five  principal  submatrices  are  2.36*  10“5, 
5.23-10-5,  5.3310-5,  5.23  10~5  and  2.36  10“5.  The  effect  of  encountering  a sequence  of  ill-conditioned  submatrices 
is  felt  later  in  the  elimination  process  and  is  manifested  by  a severe  loss  in  accuracy  in  the  subsequent  columns  of  R 
and  L.  Hence,  some  form  of  restoration  is  required  for  high  accuracy.  A way  to  bring  back  the  lost  information  is  to 
recompute  the  most  recent  columns  of  R by  solving  a Yule- Walker  problem  as  in  (1.4),  utilizing  the  decomposition 
we  already  have  at  hand.  In  this  example,  columns  10  and  11  of  R are  recomputed,  so  that  again  the  process 
restarts  from  a new  and  accurate  point.  The  results  are  presented  in  Table  5.  Notice  that  the  two  factor  matrices 
produced  by  BM  algorithm  are  nearly  singular. 
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Table  5.  Error  Behavior  for  Example  5 


BM 

Pivoted  BM 

k{L) 

k(R) 

||  x - z(0)  ||2 

II  * - *(1)  lb 
||  *-*(2)  ||, 

Il*-i(4)ll2 

#Refine. 

3.93  • 10+13 

3.93  • 10+13 

6.58  • 10“2 

2.93  • 10"5 
8.82  • 10"8 * * * 

1.58  • 10"13 
4 

1.33  • 10+3 
1.24  • 10+3 
3.29  • 10"4 
2.11  • 10-9 
4.18  • 10"14 

2 

Example  6 . We  extend  the  previous  example  to  size  50  x 50  by  appending  random  numbers.  For  this  example, 
the  Hankel  matrix  is  moderately  ill-conditioned  with  k(Hq)  = 1.97  • 10+3.  The  results  are  shown  in  Table  6. 


Table  6.  Error  Behavior  for  Example  6 


BM 

Pivoted  BM 

k(L) 

4.01  • 10+13 

1.33  ■ 10+3 

k(R) 

4.01  • 10+13 

8.71  • 10+4 

||*-*(0)  lb 

5.86- 10" 1 

2.44  ■ 10"3 

II  X - xW  lb 

2.39  • 10"3 

4.54  • 10-8 

||z-*(2)  lb 

7.46  • 10"6 

9.86  • 10-13 

||z-x(3)  lb 

j|  * - £(5)  jj2 

6.87  • 10"8 
1.17  • 10"12 

5.58 -10- 14 

# Refine . 

5 

3 

8.  Conclusion 

We  believe  that  the  Berlekamp-Massey  algorithm  works  well  when  the  Hankel  matrix  is  positive  definite  and 
well-conditioned,  so  that  none  of  its  principal  submatrices  is  ill-conditioned,  and  no  pivoting  is  necessary.  In 
general,  consider  the  Hankel  matrix  as  a moments  matrix  with  respect  to  certain  weights.  These  weights  are  not 
necessarily  all  positive,  and  thus  we  may  need  to  deal  with  Hankel  matrices  that  do  not  have  positive  definite 
property.  Strategies  such  as  pivoting,  normalization  and  gap-jumping  are  required  in  this  case. 

For  the  previous  examples,  we  adopt  a scheme  that  combines  both  pivoting  and  normalization.  To  generate  a 
new  column  of  the  triangular  factor,  say  column  i,  we  combine  it  with  columns  i — 1 and  i — 2.  Since  normalization 
is  performed  after  each  new  column  is  generated,  and  column  i is  a shifted  version  of  column  i — 1,  so  all  the  three 
columns  have  a 1 as  the  leading  nonzero  element.  Therefore,  no  pivoting  is  performed  in  the  first  phase  of  column 
combination.  After  columns  i and  i — 2 are  combined,  pivoting  takes  place  in  the  second  phase,  in  which  columns 
i and  i — 1 are  combined.  Both  steps  are  crucial  to  the  stability  of  the  procedure.  Pivoting  prevents  multipliers 
from  being  too  large,  while  normalization  keeps  the  norm  of  the  columns  from  underflowing. 

Since  we  remove  the  constraint  of  positive-definiteness,  a well-conditioned  Hankel  matrix  may  have  several 
ill-conditioned  submatrices.  The  choice  of  the  threshold  is  a subtle  issue,  and  from  the  previous  examples,  we  see 
that  a compromised  threshould  value,  such  as  10  • y/e-  ||  H ||2,  may  be  a good  choice. 

Whenever  there  is  a sequence  of  ill-conditioned  principal  submatrices,  i.e.,  a gap,  we  simply  shift  up  the 
previous  column  of  L until  the  next  well-conditioned  submatrix  is  encountered.  Thus,  we  avoid  the  computation 
within  the  gap  by  “jumping  over”  it.  Two  columns  of  L are  recomputed  right  after  the  gap,  so  that  the  errors  in 
the  factorization  caused  by  the  jumps  are  confined  within  the  gap  and  do  not  propagate  to  the  succeeding  columns. 
Therefore,  after  a few  steps  of  iterative  refinement  at  the  end,  all  the  errors  in  the  solution  (not  the  decomposition) 
will  be  corrected  and  the  solution  will  be  accurate  to  machine  precision. 
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Another  possibility  to  deal  with  the  gap  is  to  perform  a LU  decomposition  with  partial  pivoting  instead  of 
jumping  over  it.  But  the  worst  case  for  this  approach  requires  0(p2n)  operations  to  decompose  the  part  of  the 
matrix  corresponding  to  the  gap,  in  contrast  to  0(pn)  operations  for  the  gap-jumping  approach,  where  p denotes 
the  size  of  the  gap. 

Acknowledgements 

A.  W.  Bojanczyk  was  supported  in  part  by  the  Army  Research  Office  under  contract  DAAL03-90-G-0092. 
T.  J.  Lee  and  F.  T.  Luk  were  supported  in  part  by  the  Army  Research  Office  under  contract  DAAL03-90-G-0104, 
and  by  the  Joint  Services  Electronics  Program  under  contract  F49620-90-C0039. 

References 

[1]  E.  R.  Berlekamp,  Algebraic  Coding  Theory , McGraw-Hill,  New  York,  NY,  1968. 

[2]  T.  F.  Chan  and  P.  C.  Hansen,  “A  stable  Levinson  algorithm  for  general  Toeplitz  systems,”  CAM  Report  90-11, 
Computational  and  Applied  Mathematics,  University  of  California,  Los  Angeles,  CA,  1990. 

[3]  E.  Jonckheere  and  C.  Ma,  “A  simple  Hankel  interpretation  of  the  Berlekamp- Massey  algorithm,”  Linear  Alg. 
Applies vol.  125  (1989),  pp.  65-76. 

[4]  S.-Y.  Kung,  “Multivariable  and  Multidimensional  Systems:  Analysis  and  Design,”  Ph.D.  Dissertation,  Depart- 
ment of  Electrical  Engineering,  Stanford  University,  Stanford,  CA,  1977. 

[5]  N.  Levinson,  ‘The  Wiener  RMS  (root-mean-square)  error  criterion  in  filter  design  and  prediction,”  J.  Math. 
Phys.f  vol.  25  (1946),  pp.  261-278. 

[6]  J.  L.  Massey,  “Shift  register  synthesis  and  BCH  decoding,”  IEEE  lira  ns.  Inform.  Theory , vol.  IT-15  (1967), 
pp.  122-127. 

[7]  D.  R.  Sweet,  “Numerical  methods  for  Toeplitz  matrices,”  Ph.D.  Dissertation,  Department  of  Computing 
Science,  University  of  Adelaide,  Australia,  1982. 


477 


